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Abstract Classical analysis of variance requires that model terms be labeled as 
fixed or random and typically culminate by comparing variability from each batch 
(factor) to variability from errors; without a standard methodology to assess the 
magnitude of a batch's variability, to compare variability between batches, nor to 
consider the uncertainty in this assessment. In this paper we support recent work, 
placing ANOVA into a general multilevel framework, then refine this through batch 
level model specifications, and develop it further by extension to the multivari- 
ate case. Adopting a Bayesian multilevel model parametrization, with improper 
batch level prior densities, we derive a method that facilitates comparison across 
all sources of variability. Whereas classical multivariate ANOVA often utilizes a 
single covariance criterion, e.g. determinant for Wilks' lambda distribution, the 
method allows arbitrary covariance criteria to be employed. The proposed method 
also addresses computation. By introducing implicit batch level constraints, which 
yield improper priors, the full posterior is efficiently factored, thus alleviating com- 
putational demands. For a large class of models, the partitioning mitigates, or 
even obviates the need for methods such as MCMC. The method is illustrated with 
simulated examples and an application focusing on climate projections with global 
climate models. 
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21 1 Introduction 



Identifying and comparing variability among several factors is a fundamental task of sta- 
tistical analysis. From initial exploratory steps, to model testing, analysis of variance 



24 plays a vital role in the practice of statistics. Gelman (2005) has outlined a general 



25 ANOVA methodology that fits a wide range of models and summarizes results in a man- 



26 ner that facilitates interpretation across different sources of variation. Gelman and Hill 



27 (2006) elaborate further, describing the methodology in terms of a multilevel model. One 



important contribution of this approach is in providing summaries that are more con- 
structive than conclusions based on hypothesis tests. This framework is seeing usage by 



ijresearchers from diverse fields e.g. ecology (Qian and Shen, 2007), genetics (Leinonen 



et al] 2008), and climate (Sain et al. 2011). 



This paper presents a method that adopts the multilevel approach towards ANOVA, 
then extends it to multivariate settings, yielding a multilevel multivariate analysis of 
variance (MMANOVA) methodology. The strategy of initially treating all sources of 
variability in similar regard is naturally a part of the new method. We address known 



issues of applying constraints in variance analyses (Nelder, 1977, 1994, 1999) by con- 



straining batch levels so that prior distributions are improper. The extension is seen as 
a valuable contribution, since multivariate ANOVA can further obfuscate model speci- 



fication, interpretation, and computation. Kaufman and Sain (2010) offer an analysis 



of variance procedure that aids in model specification and interpretation, but requires 
computationally demanding MCMC steps. Recent methods involving approximations. 



e.g. integrated nested Laplace approximations (INLA) (Lindgren et al., 2011; Rue et al. 



2009), have allowed for MCMC to be eliminated in many cases, significantly reducing 



computational demands. While the range of problems for which such methods are appli- 
cable is wide, the focus is not typically on variance parameters. Thus, the contribution 
of our work is in promoting a general ANOVA methodology. We accomplish this by 
supporting recent work with the same goal, by refining the model specification, and by 
extending this to multivariate cases. 



2 



In Section [21 we review the ANOVA formulation of Gelman (2005) and summarize 



50 some details of the approach. We then formally extend the idea to the multivariate 

51 case, discussing technical and computational aspects. Section |3] provides a demonstration 

52 of the method through a simulation example and through a climatological application 

53 using data from future climate projections given by several atmosphere-ocean general 

54 circulation models and global emissions scenarios. In Section |4] we discuss extensions and 

55 further computational benefits that are possible when the method is applied to common 

56 high-dimensional problems. 

57 2 Analysis of Variance 

58 Analysis of variance is widely accepted as a means of partitioning variability in a manner 

59 which allows it to be attributed to various factors. An important initial step in the 

60 analysis is considering each factor of the model to be fixed or random. This step, necessary 

61 in the classical setting, raises enough issues that statisticians have been obliged to address 



62 the "mixed models controversy" (Lencina et al, 2005 Voss, 1999). One might conclude 



63 that a consensus has still not been reached, given that John Nelder deemed it necessary 

64 to reiterate the requisite points of constraints and marginality over such a long period of 



65 time, beginning with Nelder (1977) and most recently Nelder (2008). 



The rest of the section outlines a recent attempt by Gelman (2005) towards a more 



67 universal ANOVA methodology, then refines and extends the approach to a multivariate 

68 context. 



2.1 Multilevel ANOVA 



70 A fundamental contribution of the hierarchical regression approach to ANOVA employed 



71 by Gelman (2005) has been to indiscriminately consider all components in a model as 



72 random, thereby facilitating comparison across all sources of variability. The terminology 

73 is useful in supporting the indiscriminate nature of the method. The word batch is 

74 applied to all terms in the model, e.g. overall mean, factors, nested terms, interactions. 



75 etc. The nature of the variabihty from a batch is further distinguished. The distinction 

76 traditionally made by random and fixed effects is instead addressed by considering a 

77 batch's super and finite population variance. We now summarize the recent shift in 

78 ANOVA as a methodology in terms of a univariate linear model. 

79 2.1.1 Model Parametrization 



80 Following the notation of Gelman (2005 ), observations Yi,i = 1, . . . ,n are stated in terms 

81 of the additive decomposition 

B 



or the alternative regression formulation 



3 

6=0 



(b) 
b ; 
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with denoting individual batch levels and Xij denoting explanatory variables. The 



84 regression formulation could be used for the additive decomposition, with explanatory 

85 variables set to either or 1. Batch indices b = and b = B will often correspond to an 

86 overall mean, fi, and to measurement errors, e,, respectively, so that no = 1, ns = n. An 

87 individual batch is referenced by 6 = 0, . . . , -B and consists of rii, levels. Individual levels 

88 of a batch are denoted as /3[^\ . . . , /Sn,^ with jj, . . . , replicating the levels so that each 

89 observation is associated with exactly one batch level. We acknowledge that the additional 

90 level of subscripts and superscripts may seem contrived to some, although it is necessary 

91 for the general case. In practice the number of batches in the model is reasonable so 

92 that this can be avoided, as done in Section [3j Additional sub or superscripts can often 

93 be dropped. For example, a batch in its entirety is denoted by /3^''^ = . . . ,/3n^}- 

94 Given a batch, and an assumed distribution on model errors, a conventional fixed effects 

95 analysis often corresponds to the test Hq : Pj'^^ = for j = 1, . . . ,nf,. While for a random 

96 effects model, assuming the nj, levels of each batch b to be modeled as Gaussian 

(3fU al NiO, al), j = l,...,n,, (3) 



97 a test for significant batch variation would be Hq: = 0. Alternatively, the proposed 

98 methodology identifies two representations of variation of a given batch. The superpopu- 

99 lation variance, a^, corresponds to the variance of all potential, possibly infinitely many, 

100 levels of a batch. The finite-population variance represents variability of the specific set of 

101 batch levels that have been realized. Super and finite-population variances can be roughly 

102 related to the random effect variance component estimate, and the fixed effect within- 

103 group sum of squares, respectively. As an example, consider batch b and its vector of batch 

104 levels, f3^^^ = ■ ■ ■ .Pn})'^ with Cb constraints. Then the degrees of freedom are z/^ = 

105 rib — Cb, and the finite-population variance is sf = (in,, — {CbC'[)'^Cb) /3^^\ 

106 where I^^ is the x Ub identity and is the Cb x n?, constraint matrix such that 

107 Cb(3^''^ = 0. Variance component estimation is made by decomposing the variance of the 

108 batch level estimates, H = Var(/3('')) = Var{E(^('') | + E{Var(/3('') | /?(''))}, into the 

109 sum of the superpopulation variation, a^, plus the variability of the batch level estima- 

110 tions, Vb;estimation- The choscu estimate of the superpopulation variance, in this case the 

111 met ho d-of- moments estimator, is then a"^ = Vb — Vb:estimation, where, H = ^ Y^^LiiPf^Y^ 

112 and Vb:estimation = X]fce/(6) n^^k includcs supcrpopulatiou variances from other batches 

113 that enter into variability of the batch level estimates, indicated by the set I{b). At a 

114 minimum, Vbiestimation will include a^, the estimated error variance. For a large class of 

115 multilevel, hierarchical models, this strategy allows for all terms to be treated as random, 

116 and for their variabilities to be assessed. 



117 2.1.2 Confirmatory Procedures 

118 In regards to more inferential procedures, either a frequentist or Bayesian direction can 

119 be taken. In the frequentist case, an inverse-chi-square distribution, is employed 

120 to assess uncertainty in the superpopulation variance a^; since -^sl/a^ is chi-square dis- 

121 tributed with z/f, degrees of freedom. For batches with X]fcG/(f)) nt^^' i^icluding more than 

122 only the error variance a^, then a linear combination of inverse-chi-square distributions. 



123 J2i ^iXu^y is required, as described at the end of the previous section. As Gelman (2005 ) 
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124 states, these linear combinations may be dealt with directly, although simulation is often 

125 more straightforward. The simulation, described therein, is carried out by: 1) Obtain 

126 R simulated raw variances, Vb for each of the B batches in the model with a random 

127 variable that is proportional to and corresponding degrees of freedom; 2) Calculate 

128 superpopulation variances using = max(0, H — Vbiestimation); 3) Simulate batch levels 

129 using newly generated superpopulation variances; 4) Calculate sample variances of each 

130 batch. This procedure then yields a (posterior) sample of superpopulation variances, 

131 batch levels, and finite-population variances, corresponding to the final three steps. 

132 A strict Bayesian approach requires additional prior specifications, but yields posterior 

133 distributions of the superpopulation variances. However, this distinction, between the 



134 two schools of thought, can be seen as purely semantic. As Gelman (2005) states, "given 

135 cr^, the parameters /Sj''^ have a multivariate normal distribution (in Bayesian terms, a 

136 conditional posterior distribution; in classical terms, a predictive distribution)". Thus, 

137 assessing uncertainty in the finite-population variances, sf, is the same. In both cases, 

138 either batch levels themselves are simulated, or the distribution of is approximated 

139 with an appropriate chi-square random variable. 

140 It should be clear that the uncertainty surrounding superpopulation variance param- 

141 eters will typically be greater than that for finite-population variances. Intuitively, this 

142 is because superpopulation variances describe variability of levels that have not yet been 

143 realized. 

144 2.2 Multilevel Multivariate ANOVA 

145 Typical multivariate analysis of variance strategies rely on the distribution of the de- 



146 terminant of sums of squares matrices, i.e. Wilks' lambda distribution (Mardia et al. 



147 1979, p. 335), and culminate in value related conclusions. Hence, it does not easily 

148 facilitate inclusion of random effects, and thus no comparison across these effects. Using 

149 a Bayesian approach, we now derive a general multivariate methodology that seeks to 



150 provide results similar to those of Section 2.1.1 The method partitions variability by 
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151 batch in an efficient manner and further factors the posterior by batch into a batch's 

152 superpopulation covariance posterior, and batch levels posteriors. In addition, we handle 



153 the issue of constraints in a way that does not commit what Nelder (1994) has called 

154 one of the false steps of linear models. Rather, the constraints are implicit, yielding im- 

155 proper batch level prior distributions. Another point which must be mentioned is that 

156 of matrix parameter estimation. Although the common limitations of covariance matrix 

157 estimation and modeling are issues that must be dealt with, it is important to ffist focus 

158 on, and refine multivariate ANOVA for familiar settings. In Section |4] we discuss issues, 

159 modifications, and implications of the method when higher dimensional data is used. 

160 2.2.1 Multivariate Model Parametrization 

161 Consider d- dimensional multivariate observations such that batch levels are vectors and 

162 batch variances are covariance matrices. Namely, (jT])-(|3]) now contain vectors Yj, /3^b\ 

163 matrices Xj j, and covariance matrices are S^, all of appropriate dimension. We adopt the 

164 strategy from the previous section, in indiscriminately considering all terms as potentially 

165 possessing variability. The multivariate analogue of ([!]) with Gaussian errors is 

Y,\{(3^P},,l:,r^^^,(^^f3^^, j , ^ = l,...,n, (4) 

166 and the remainder of the Bayesian model specification is given by 

/3f |S5~Ar,(/3f), S,), J = 1,..., rife, (5) 

Vb = ^b+-^e\^e--W-\<ffb, f^b), (6) 
n 

~ W-\^,, (7) 

167 for 6 = 0, . . . , -B — 1 and with the inverse-Wishart distribution denoted by W~^. Batch 

168 indices 6 = and b = B respectively correspond to the intercept, or overall mean term, 

169 n, and to measurement errors, e^. For notational convenience we will refer to and 

170 rather than Sq and S^. Typically zero-mean batch level priors are assumed; that is, 

171 f3Q^ = 0. Setting = 0, ki, = yields the noninformative prior p{\Jb) |Ub|^'^'^+"'^^/^. 



172 Because inverse- Wishart support is given by the set of all positive definite matrices, ^ is 

173 referred to as a constrained inverse- Wishart distribution, since Uf, — is required to be 

174 positive definite. This covariance parametrization, U;, = Sf, + — S^, has previously been 



175 utilized in the context of multivariate random effects by Everson and Morris (2000), for 



176 which they develop efficient methods of sampling. Also, when only error terms contribute 



177 to the variance of batch level estimates, \Jb is analogous to Vb of Section 2.1.1 with 

178 I{b) = {B}. 

179 The choice of inverse- Wishart covariance priors, ^ and ([7|, has been made to balance 

180 the complexity of model specification with implementation and computation. However, 

181 given the considerable amount of research of covariance priors, there are other options 



182 available. Daniels (1999) and Daniels and Kass (2001) examine covariance priors that 



183 emphasize uniform shrinkage of the eigenvalues. Although informative covariance priors 

184 may be necessary in many cases, the usage of such priors will impact the computational 

185 demands required. 

186 2.2.2 Posterior Distributions 

187 Without additional specification, (|4])-([7]) yield an inadequate posterior. In an MCMC 

188 setting this may manifest itself by failure to converge, due to drifting in the parameter 

189 space. From a classical point of view, estimating the set of all batch levels, {/3jt^} would 

190 require additional constraints. The inclusion of similar constraints in the Bayesian model 

191 allows for a closed form of the posterior, as well as for factorization between batches. 

192 Degrees of freedom for each batch are then accounted for in the corresponding batch 

193 covariance posterior. This parametrization is also beneficial in terms of computation 

194 since batches are conditionally independent of one another. Using a vectorized form of 



195 the model, Y = (Yf , . . . , Y^)^ G M"*^ and /S'^^^ = {/Sf^^, /S^^^'^f E M""'^, is convenient 

196 for the development. The constraint Cb/?*-^-* = G M.'^'''^, where there are Cb constraints, 

197 combined with (|4]), ([s]), is now 
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Y I {(3^%, Ctf3^'^ = 0, ~ Af, (j^f^^'^^ I- ® ' (8) 
/J^") I Cb/3('') = 0, Sfc ~ A/;,,(l„, ® fife), (9) 

198 with denoting the Kronecker product. The rank-deficient Qb, due to the constraint, 

199 causes Q to be improper. To derive this improper distribution begin with the uncon- 

200 strained and vectorized form of ([s]), which has covariance fij, = (g) S;,. The density is 

201 then stated through a decomposition of the precision, = FAF"^ ® S^^^ Assuming Cb 

202 constraints, the rank deficiency is introduced by removing the corresponding number of 

203 eigenvalues from the diagonal matrix A, e.g. A = diag(0, . . . , 0, Xc^+i, ■ ■ ■ , AnJ, leading to 

204 Qb = FAF"^ (g) S^^. This method of addressing linear constraints is useful for other gen- 

205 eral improper distributions and intrinsic Gaussian Markov random fields, as illustrated 



206 by Rue and Held (2005). Let |Q|* be the pseudo determinant of a singular matrix, that 

207 is, the product of its non-zero eigenvalues. In the case of the identity matrix being used 

208 as the first matrix term in the Kronecker product, the eigenvalues are one, thus densities 

209 of the likelihood, ([s]), and of batch level priors, ^ are 

/ n 

p{Y I {/3f }fe, S,) (x|S,|-"/2exp 5^(Y, - Y,fi::\Y, - Y,) 

\ i=i 

^ B—l rii, 

6=0 j=l ^ 

pifB'-"^ I S,) oclQfeiy^exp - 1 ® (3^^Y Qt{(3'^''^ - 1 ® /^i'^)) 

oc|E,|-(---)/^exp (-1 f; -/3f^)^S-(/3f -/3(^))) , 

210 where ^ denotes least-squares estimates. For orthogonal batches, the full posterior from 

211 these densities and from batch covariance prior densities, can be conveniently factored 

212 into 



B-l 



p(S„ {Sfe, /3W}f^-i I Y) = p(S, I Y) n (3^'^ I Y, S,). (10) 



b=0 



213 Each joint batch density, f3^'^^ \ Y, S^), is then factored further. Using known matrix 

214 identities, e.g. (A^^ + B"^)"^ = A(A + B)"^B, the identity 

{x - sfV-\x -s) + {x- tfY-\x - t) 

(11) 

= tT {(U + y)'\s - t){s - tf) + {x -P-^mfP{x -P-^m), 

215 is derived, where P = U^^ + V^^, m = U^^s + V^^i. Accounting for model constraints. 



216 together with (11), the batch superpopulation posterior and the batch levels are found 



217 through the decomposition of quadratic forms of batch levels and least squares estimates 

Tli, 
j = Cb + l 

218 where P^ = Sr^ + ^Sr\ mf^ = ^S-^3f^ + Sr^/3^\ and tr(-) denoting the trace 

219 operator. Additionally, B;, = X]j=c6+i(/^i ^ ~ (^o^)iP\ ~ f^'o^Y ■> is analogous to a matrix 

220 sums of squares of the unconstrained batch level estimates that has been adjusted by the 

221 prior mean. The full joint posterior is then factored as 

B-l 

{S,, /3W}fro^ I Y) = I Y) J] I Y, p{(3^'^ \ Y, E„ S,), (12) 

b=0 

222 where the product denotes batch posterior independence, and thus no need for computa- 



223 tionally intensive MCMC procedures. The corresponding distributions of (12) are 

(n B-l \ 

+ J](Y,-Y,)(Yi-Y,)^, K, + n-Y,nA, (13) 
i=l fe=0 / 

Sfc + ^S, I Y, ~ W-^ (*6 + Bfe, K, + n,- c) , (14) 
n 

/3f I Y,S„S,~ <^ V ^ - ; (15) 

(P,-'mf , P-I) J = c, + 1, . . . ,n,. 

224 Batch levels (15), which reflect both free and constrained parameter estimates, can then 

225 be sampled and adjusted accordingly to obtain posteriors for finite-population covari- 

226 ances, S^. Recall finite-population parameters focus on observed levels of a batch, not on 

227 all potential unobserved batch levels. 
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2.2.3 Covariance Posteriors 



229 In all cases thus far covariances are assumed to be of full rank. Hence, improper 

230 covariance posteriors will be due only to an insufficient number of observed levels of the 



231 batch. Diaz-Garcia et al. (1997) offer a comprehensive look at all possible cases of im- 



proper Wishart distributions and following their terminology this would be classified as 



233 a pseudo- inverse- Wishart. Uhlig (1994) as well as Srivastava (2003) consider sampling 



234 with a pseudo-singular- Wishart distribution. Extending their work to inverse- Wishart 

235 distributions is one method for dealing with moderate discrepancies in the number of 

236 observed levels, ri;, < d. Addressing cases in which rib <^ dis discussed in Section |4] Even 

237 in the remaining case, > d, simulation from a posterior is not always efficient. Because 

238 support of the inverse- Wishart posterior requires positive definiteness in two respects. 



Hp > 0, and S 



/3 



rib 



> 0, the usual method of rejection sampling from an inverse- 



Wishart distribution is not always practical. Everson and Morris (2000) describe a more 



241 computationally efficient method to maintain positive definiteness through a Cholesky 

242 decomposition and maintaining positive eigenvalues while the sample realization is gen- 

243 erated. 



244 2.2.4 Analysis Results 

245 Multivariate sources of variability do not always yield a single, clear criterion that in- 

246 dicates the greatest contributor to overall variability. For scalar variance components, 

247 > al is clearly interpreted, however, due to the partial ordering of positive definite 

248 matrices, the analogous statement on covariance matrices is not useful. In other words, 

249 there is not a single, obvious comparison that can be made to determine which of two 

250 covariances are "greater" . Depending on the setting, there may exist an adequate scalar 

251 that sufficiently summarizes covariance characteristics. For the volume of ellipsoidal con- 

252 tours the determinant, |S|, achieves this, while in other cases the sum of all entries, 

253 l^Sl, or the sum of the marginal variances, tr(I]), may be appropriate. Scalar criteria 

254 with corresponding uncertainty intervals then allow multivariate sources of variability to 
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255 be directly compared. Mardia et al. (1979) employed the determinant and trace, which 

256 correspond to the product and sum of eigenvalues; referring to them as the generalized 

257 variance and total variance, respectively. We have amended the terminology, since inclu- 

258 sion of the sum of all matrix elements requires further delineation, by denoting 1-^Sl as 

259 the total variance, and tr(S) as the total marginal variance. 

260 Effectively relaying results of analysis of variance is one of the motivating factors 



261 that Gelman (2005) cites. The classical table of values does not yield any indication 

262 of batches with the largest variances, nor are the required assumptions on other batches 



263 clear. Nelder ( 1999 ) discusses many of the issues of over-reliance upon the p- value, and its 

264 ineffectiveness as a tool for communicating results. Uncertainty intervals are the default 

265 choice for presenting results. Visual plots are convenient since they facilitate simultaneous 

266 comparison of the relative variability contributions, their magnitudes, and the magnitude 

267 of the uncertainty in the estimates. For direct comparison of batches h and h\ statements 

268 of the form P((7(Sb) > 5'(Sfc/) | Y), utilizing an arbitrary matrix criterion (/(■), may also 

269 be found. 



270 3 Examples 

271 This section covers two examples of the outlined methodology to carry out confirmatory 

272 procedures on multivariate data. The first is a toy example in which output and sum- 

273 maries are given in order to provide further insight. The second example utilizes global 

274 averages of temperature and precipitation predictions using 13 atmosphere-ocean gen- 

275 eral circulation models and 3 global greenhouse gas emissions scenarios that have been 

276 identified by the Intergovernmental Panel on Climate Change. 



278 



279 



3.1 Simulation 

The process Yjj = ^ + CKj + e^j G W^, i.e. -B = 2, (i = 3, is used to illustrate the method 



of Section 2.2 To generate the three-dimensional observations, parameters Sq,, S^, and 



280 /X are fixed. An individual simulation is then performed by generating CKj and ejj,i 
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281 1, . . . , no,, J = 1, . . . , n^, from mean-zero multivariate Gaussian distributions with their 

282 respective fixed covariances. The mean term fi is added to generated data resulting in 

283 n = UaTi^ observations. 



284 Using (13)-(15), posterior distributions for these covariance criteria are obtained for 

285 three distinct simulation scenarios. The three scenarios can be explained by the am- 

286 biguous description that variability introduced by batch ck is greater than (case 1), less 

287 than (case 2), or comparable to (case 3) variability introduced by error batch e. More 

288 specifically, covariance matrices are decomposed into a vector of marginal standard devi- 

289 ations s and a correlation matrix R, e.g. Sq, = diag(sa)RQ,diag(sQ,). For all simulation 

290 correlation matrices are fixed. Batch levels CKj have the unique correlation structure 

291 (Rq,)i,2 = 0.3, (Rq)i,3 = 0.1, (Ra)2,3 = 0.5. Errors Sij have the correlation matrix, R^, 

292 with autocorrelation structure (Re)i,j = p'*"''', and p = 0.2. Error marginal variances are 

293 additionally held constant at 1 over all simulations, = R^, so that only is distinct 

294 for each case. 

295 The objective of the analysis is to assess the relative variability introduced by batch 

296 a and batch e, as well as the uncertainty in the assessment. Further, this is to be done in 

297 an appropriate multivariate context. Figure [T] displays the results of two simulation runs 

298 under the first scenario (case 1) using the determinant. In one simulation, the number 

299 of batch level realizations are ria = 5,n^ = 3 and in the second Ua = 8,n^ = 5, which is 

300 to say less vs. more data. The left-most graph of Figure [T] displays uncertainty intervals, 

301 with narrow lines denoting 95% quantiles, thicker lines 50%, and a vertical tick placed 

302 at the median. The upper set of intervals, which are intuitively wider, correspond to less 

303 data, while the lower set of intervals correspond to more. By vertical comparison of the 

304 uncertainty intervals, we see that for = 5, = 3 all intervals overlap, and hence no 

305 distinction can be made between the sources of variability. For ria = 8,n^ = 5 however, 

306 there is no overlap of the uncertainty intervals, suggesting that both superpopulation 

307 variability, and the finite-population variability are greater than error variability. Ad- 

308 ditionally. Figure [T] offers a diagnostic look at covariance posteriors. The center figure 
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309 displays ellipses from first two principal components corresponding to 2.5%, 50%, and 

310 97.5% determinant-ordered percentiles of the posterior distributions, specific values of 

311 which can be seen on endpoints of the uncertainty intervals. The right-most graph of 

312 Figure [T| which shows all 1000 ellipses from each posterior distribution, offers a look at 

313 the size, shape, and orientation, for an overall comparison of the uncertainty in the batch 

314 covariances. Size renders an idea of the magnitude of the marginal variances. The dis- 

315 parity in size between the ellipses corresponding to Sq, and suggest that the marginal 

316 variances of the former are greater than those of the latter. Through shape, one may 

317 glean some insight into batch covariance dependence. Lastly, the orientation, or varying 

318 orientation, suggests the uncertainty of the dependence, e.g. as the orientation of the 

319 ellipses corresponding to fluctuate greatly, there is not much that can be said about 

320 its dependence structure. 

321 For comparison, consider a frequentist approach to a simplified form of the problem. 

322 Beginning with Xj ~ A/'d(0, S), z = 1, . . . , ra, it is known that U = ^^(Xj— X)(X,;— X)"^ ~ 

323 PVrf(S,n — 1), from which we derive the distributions 

d 

\v\^m-l[xi-i-d+i, 



i=l 



l^Ul~l^Sl-x^_i, 

d 

tr(U) ~ 5^ A, ■ xl-i, Ai, . . . , A, = eig(E). 

1=1 

324 The first two offer pivots and thus allow for closed form expressions that yield confidence 

325 intervals for values of interest |S|, 1-^Sl. The confidence interval for tr(S) is based on 

326 the normal approximation that matches the first two moments, E{tr(U)} = (n— l)tr(S), 

327 and Var{tr(U)} = 2{n — l)tr(S^). This approximation has been chosen in the spirit of 



328 moment matching approximations used by Imhof (1961) as applied to quadratic forms 

329 of random vectors. Note that these confidence intervals assume that the CKj are directly 

330 observed, which is not the case for our proposed method. Rather, the classical methods 

331 shown are included only for comparison. 

332 To gain insight into the coverage success and uncertainty interval widths, we have 
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333 carried out S = 100 simulations over different values of and each of the different 

334 scenarios of variability sources. For all simulations the number of replicates at each level 

335 is fixed at = 15. Results in all 3 cases were relatively similar with respect to coverage, 

336 noting however that uncertainty interval widths increase as the magnitude of variability 

337 increases. Thus, only the scenario in which variability sources are comparable (case 3) 

338 has been shown (Figure [2]). Coverage and interval widths for Sq, and Spreq should be 

339 compared as they correspond to the same true, unknown covariance. Despite the fact 

340 that the methodology does not assume realizations of a to be observed directly, but 

341 rather indirectly through Y, results are comparable to the frequentist approach outlined. 




' — I 1 ] ] ] ' o ^1 1 — — I ^1 1 ^ 1 

10"^ 10" 10^ 10^ 10^ 10 20 30 40 10 20 30 40 

Figure 1: Simulation for two distinct sample size pairs, Ua = 5,n^ = 3, and Ua = 8, = 5. 
In both, batch a varies greater than the errors, with = (v^, V^, V^)"^- Uncertainty 
intervals with thin narrow lines denoting 95% uncertainty level endpoints, thicker middle 
portion denoting 50% endpoints, and vertical tick at median. Larger red vertical ticks 
denote true parameter values. Higher-positioned, wider intervals correspond to smaller 
batch sample sizes, while lower-positioned, narrower intervals correspond to larger batch 
sample sizes (left). For second set of sample size pairs, ellipses from first two principal 
components are displayed for 2.5%, 50%, and 97.5% determinant-ordered percentiles of 
posterior distributions when 1000 posterior realizations have been drawn (center). All 
1000 posterior ellipses (right). 
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Figure 2: For Ua = 5, . . . , 100, coverage is estimated with S = 100 simulations in which 
the sources of variabihty are comparable (case 3). Coverage estimates for S^, S^, and 
where grey line denotes 95% nominal coverage (top). Average uncertainty interval 
widths on log scale (bottom). 

3.2 Application 

In this example our methodology is applied to a bivariate dataset of global temperature 
(Celsius) and precipitation (mm/day) for 9 decadal averages of boreal summer months, 
June, July, August, during the remaining century. The first batch in the model con- 
sists of 13 levels, each representing a single atmosphere-ocean general circulation model 
(AOGCM) developed by several international climate research institutions as part of the 



348 CM1P3 project (Meehl et a/., 2000) in the framework of the Fourth Assessment Report 



(AR4) for the Intergovernmental Panel on Climate Change (IPCC). The second batch 
covers 3 greenhouse gas emissions scenarios that have been defined by the Special Report 



ipn Emissions Scenarios (SRES), which are identified as AlB, A2, and Bl (Nakicenovic 



and Swart, 2000). One fundamental objective of the analysis is then to compare how 



these factors contribute to overall variability of global climate averages, how they relate 
to one another, and what the uncertainty of this assessment is. 

Bias and dependence among climate models is an issue that has more recently begun to 



Jbe examined further, beginning with Tebaldi and Knutti (2007), Jun et al. (2008), Knutti 



et al. (2010), and references therein. Despite this, we adopt the statistical assumption 
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358 that has traditionally been used when with working with sets of AOGCMs, which is 

359 to assume that they are independently drawn from a common process representative of 

360 true climate characteristics. Using this assumption our approach can be seen as a useful 

361 exploratory tool, and may be further adopted to address contrasts of batch levels, and 

362 thus to identify similar batch levels. Preliminary analysis steps have suggested the model 



ijt 



fJ'Q + "O.i + /3oj + 7ij + Xi^tfJ'i + Xi^t0^l,i + Xi^tf^^j + X2,tfJ'2 + f^ijt, (16) 



363 where ? = 1, . . . , tIq, = 13, j = 1, . . . , n^^ = 3, t = 1, . . . , = 9, n = rtanpnt, and d = 2. 

364 Time covariate xi is centered such that xi^t = —4, . . . , 4, and X2 is transformed to be 

365 orthogonal to other predictors in the model. Batches of interest are AOGCM, a, and 

366 SRES, f3, and their interaction, 7. The first two are further specified as a constant effect, 

367 (Xq, fBf), as well as with respect to time, ai,f3i. 



368 Posterior distributions of batches aQ,f3Q, and 7 are derived from (14) and (15). 

369 Batches cki and f3i differ slightly as they correspond to the regression model formulation. 

370 Multivariate batch levels associated with a covariate would, in general, be multiplied by 

371 a matrix, e.g. Xi t. Using a matrix covariate, superpopulation and batch level posteriors 

372 of batch CKi are then 

/ ria \ 
Sai + I Y, ~ Yl -CaA , (17) 



oti,i I Y, Se, Sq,^ 



Nd {P^lmo,^.i, V^l) i = c„i + 1, . . . , no 

373 where = ^^{Y:tU^lt^:^^i,tr\ = + V'S and m„,,, = V-iSi,,. For 



''ai 



374 model (16) the covariate matrix is Xi_t = diag(a;i,t), and thus Vq^ = {rip X]r=i ^1 1) ""^^e- 

375 The posterior of batch is found similarly. 

376 Figure |3] suggests that AOGCM is the most distinguishing feature. Figures |4] and 

377 [5] confirm this assessment since ckq is seen as the most significant source of variability 

378 among all batches. Comparison of Figures 111 |5] also illustrate the additional uncertainty 



379 of superpopulation parameters over their finite-population counterparts. Superpopulation 
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380 covariance criteria uncertainty intervals are wide because they account for uncertainty 

381 in unobserved batch levels, particularly in the case when a small number of batch lev- 

382 els have been observed. Finite-population covariance uncertainty intervals are generally 

383 smaller, because they are concerned with variability of only the batch levels that have 

384 been realized. 






2010 2030 2050 2070 2090 2010 2030 2050 2070 2090 2010 2030 2050 2070 2090 

Figure 3: Global temperature (top) in degrees Celsius and precipitation (bottom) as 
mm/day over rit = 9 decadal intervals for Ua = 13 AOGCMs and rip = 3 emissions 
scenarios. 
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Figure 4: Superpopulation covariance uncertainty intervals using determinant, total vari- 
ance, and total marginal variance criteria shown on a log scale. Nominal coverages of 
0.95 (thin lines) and 0.50 (thick lines), and the median (vertical line) are denoted using 
quantiles of the corresponding posterior distributions. 



385 Posterior predictive distributions, often used to perform model checking and diagnos- 
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Figure 5: Finite-population covariance uncertainty intervals using determinant, total 
variance, and total marginal variance criteria shown on a log scale. Nominal coverages of 
0.95 (thin lines) and 0.50 (thick lines), and the median (vertical line) are denoted using 
quantiles of the corresponding posterior distributions. 

tics, can also be utilized to identify distinct sources of variability. The posterior predictive 
is conditional on observations with levels from each batch assumed to be, a) the same as 
those batch levels that have been observed data, b) unobserved/novel batch level realiza- 
tions. Figure [g] examines posterior density p(Yjj„^ — Yjji|{Yjj(}), the difference between 
posterior predictive distributions at the final, t = rit, and initial, t = 1, decades. Thus, 
the focus is on temporal batches, cxu', fBy,. Indices signify new, unobserved batch 
levels. Linear and quadratic terms, i-ii, I-I2 are included, although additional variability 
from these terms has been disregarded. The left-most panel, in which every batch con- 
tributes a new batch level realization, shows a large degree of variation. The center panel 
assumes that the AOGCM observed in the original data is to be used, thus variability 
from these specific batch level posteriors is included. For SRES a novel batch level is 
assumed, thus a realization utilizing the SRES superpopulation posteriors is included. 
A subset of three observed AOGCM levels has been displayed, selected so as to best 
represent the range and relative distances of their peaks. However, the high degree of 
variation introduced by the new emissions scenario level makes even these distributions 
nearly indistinguishable. In the right-most plot, using observed emissions batch levels, 
rijs = 3, the additional variability introduced comes primarily from the new AOGCM 
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403 batch level to be observed. The posterior predictive plots are particularly useful for de- 

404 termining whether the variability from each batch is due to the magnitude of the batch 

405 variability itself, or due to uncertainty in the assessment itself. 
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Figure 6: Posterior predictive distribution, p(Yi'j/nt — ^i'j'ilV^ijt}), where both 
represent new, unobserved AOGCM and SRES batch levels (left). Posterior predictive 
pO^ij'nt ~^ij'i\V^ijt}) with observed AOGCM batch levels, i — 3,9,12 (solid, dashed, 
dotted) and unobserved SRES batch level I3i j, (center). Posterior predictive p(Yj/j„j — 
Yi,ji\{Yijt}) with unobserved AOGCM batch level ai,i/ and all SRES levels, j = 1,2,3 
(solid, dashed, dotted), that have been observed (right). Density contours correspond to 
quantiles 0.05,0.25,0.75. Horizontal and vertical axes denote precipitation in mm/day 
and temperature in degrees Celsius, respectively. 



406 4 Discussion 

407 The first contribution of this paper has been in extending recent philosophical shifts in 

408 the treatment of analysis of variance to multivariate settings. New analysis of variance 

409 approaches allow appropriate parameters, e.g. super or finite population, to be used to 

410 answer the correct research question, while at the same time providing coherent model 

411 definition, implementation, and interpretation. This same fiexibility has been extended 

412 to multivariate cases; in that the researcher can guide covariance criteria choices, rather 

413 than the method determining the criterion. The second contribution has been in provid- 

414 ing a foundation for computational efficiency, which is necessary for dimension scalability. 
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415 Using improper batch level priors we have shown that it is possible to minimize depen- 

416 dencies between batch covariances. In many cases this reduces, or eliminates, the need 

417 for complex and computationally demanding analyses. 

418 Further extensions to the methodology must explicitly address increasing dimension- 

419 ality. For moderately sized dimensions d, relative to number of observations, improper 

420 inverse- Wishart distributions, and/or priors that impose particular dependence struc- 

421 tures, are possible options. For cases in which d is very large, stricter covariance as- 

422 sumptions may be employed. In the spatial context, properties such as stationarity allow 

423 covariance parameter space to be reduced, e.g. range, sill, and nugget in a spatial covari- 

424 ance function. Because simultaneous estimation of such parameters is nontrivial, some 

425 parameters are often assumed, or estimated empirically in earlier analysis steps, as in 



Furrer et al. (2007). In other cases, so as to maintain computational feasibility, spar- 



427 sity restrictions are placed on covariances (Cressie and Johannesson, 2008 Furrer et al. 



428 2006 Stein, 2008). For many such scenarios a covariance is decomposed into a correla- 

429 tion matrix and a scalar variance parameter. Our method is then carried out with the 

430 inverse- Wishart posterior density transformed through a spectral decomposition of the 

431 correlation matrix, thus allowing for efficient posterior sampling for cases in which d^ n. 

432 This extension offers an alternative to geostatistical model analyses that have previously 

433 relied on computationally intensive MCMC methods, and is the focus of current research. 

434 Other difficulties encountered are unbalanced designs and linearly dependent predictors. 

435 MCMC may be utilized for sets of dependent batch levels. Development for these cases 

436 is another area of current research. 
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